# !pip install yahooquery pandas_datareader arch
# !pip install keras-tunerPredicting stock prices with volatility and Greeks
This project pulls together a few ideas I care about: how volatility persists in markets, and whether a model that understands Black–Scholes pricing can get a bit closer to tomorrow’s close. I build a feature set around MSFT that mixes the usual technical indicators — moving averages, log returns, realized volatility — with a GARCH(1,1) estimate of volatility and the option Greeks (delta, gamma, theta) derived from it. The target is the next day’s closing price, and an LSTM is asked to learn how those features unfold over time.
The pipeline is standard but touches a lot of moving parts: fetch the data, engineer features, turn the series into look-back sequences, tune the architecture with Keras Tuner, then evaluate on a held-out test stretch. The options data gets pulled in too; it’s not used directly in the features yet, but keeping it around leaves room to learn from implied volatility later.
import datetime as dt
import pandas as pd
import yfinance as yf
import numpy as np
from pandas_datareader import data as pdr
from scipy.stats import norm
from arch import arch_model
import warnings
warnings.filterwarnings("ignore", category=FutureWarning)
# --- Data and Feature Engineering ---
def build_black_scholes_dataset(ticker_symbol, start_date, end_date, option_type='call', default_T=30/365):
tk = yf.Ticker(ticker_symbol)
data = tk.history(start=start_date, end=end_date)
data = data.reset_index().set_index('Date')
data.index = data.index.tz_localize(None) # Ensure tz-naive index
data['Log_Returns'] = np.log(data['Close'] / data['Close'].shift(1))
data['Realized_Volatility'] = data['Log_Returns'].rolling(20).std() * np.sqrt(252)
data['SMA_20'] = data['Close'].rolling(20).mean()
data['SMA_50'] = data['Close'].rolling(50).mean()
returns = data['Log_Returns'].dropna() * 100
am = arch_model(returns, vol='Garch', p=1, q=1, dist='normal')
res = am.fit(disp='off')
data['GARCH_Volatility'] = res.conditional_volatility.reindex(data.index) / 100 * np.sqrt(252)
treasury = pdr.get_data_fred('DGS10', start_date, end_date) / 100
treasury.index = treasury.index.tz_localize(None)
treasury = treasury.reindex(data.index, method='ffill')
expirations = tk.options
options = []
for expiry in expirations:
try:
chain = tk.option_chain(expiry)
opt_df = chain.calls if option_type.lower()=='call' else chain.puts
opt_df['expiration'] = pd.to_datetime(expiry)
options.append(opt_df)
except:
continue
df_options = pd.concat(options) if options else pd.DataFrame()
features = pd.DataFrame({
'Underlying_Close': data['Close'],
'Log_Returns': data['Log_Returns'],
'SMA_20': data['SMA_20'],
'SMA_50': data['SMA_50'],
'Realized_Volatility': data['Realized_Volatility'],
'GARCH_Volatility': data['GARCH_Volatility'],
'Risk_Free_Rate': treasury['DGS10'],
'Time_to_Maturity': default_T
}, index=data.index)
def bs_greeks(row):
S = row['Underlying_Close']
K = S # ATM approximation
T = row['Time_to_Maturity']
r = row['Risk_Free_Rate']
sigma = row['GARCH_Volatility']
if T <= 0 or sigma <= 0:
return pd.Series([np.nan, np.nan, np.nan], index=['Delta', 'Gamma', 'Theta'])
d1 = (np.log(S/K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
if option_type.lower() == 'call':
delta = norm.cdf(d1)
theta = (-(S * sigma * norm.pdf(d1))/(2 * np.sqrt(T)) - r * K * np.exp(-r * T) * norm.cdf(d2)) / 365
else:
delta = norm.cdf(d1) - 1
theta = (-(S * sigma * norm.pdf(d1))/(2 * np.sqrt(T)) + r * K * np.exp(-r * T) * norm.cdf(-d2)) / 365
gamma = norm.pdf(d1) / (S * sigma * np.sqrt(T))
return pd.Series([delta, gamma, theta], index=['Delta', 'Gamma', 'Theta'])
greeks = features.apply(bs_greeks, axis=1)
features = pd.concat([features, greeks], axis=1)
return features, df_options
# --- Data Preparation for LSTM with Dates ---
def create_sequences(X, y, dates, lookback=10):
X_seq, y_seq, date_seq = [], [], []
for i in range(len(X) - lookback):
X_seq.append(X[i:i+lookback])
y_seq.append(y[i+lookback])
date_seq.append(dates[i+lookback])
return np.array(X_seq), np.array(y_seq), np.array(date_seq)
def prepare_data_lstm(data, target='Underlying_Close', test_ratio=0.2, lookback=10):
data = data.dropna().copy()
dates = data.index # Preserve the date index
data = data.select_dtypes(include=[np.number])
data['Target'] = data[target].shift(-1)
data = data.dropna()
dates = data.index # Updated dates after dropping NA
X = data.drop(columns=[target, 'Target']).values
y = data['Target'].values
from sklearn.preprocessing import MinMaxScaler
scaler_X = MinMaxScaler()
scaler_y = MinMaxScaler()
X_scaled = scaler_X.fit_transform(X)
y_scaled = scaler_y.fit_transform(y.reshape(-1, 1))
X_seq, y_seq, date_seq = create_sequences(X_scaled, y_scaled, dates, lookback=lookback)
split = int(len(X_seq) * (1 - test_ratio))
return X_seq[:split], y_seq[:split], date_seq[:split], X_seq[split:], y_seq[split:], date_seq[split:], scaler_X, scaler_y
# --- LSTM Model Building with Keras Tuner ---
import tensorflow as tf
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense, Dropout
import keras_tuner as kt
from sklearn.metrics import mean_squared_error
def build_lstm_model(hp, input_shape):
model = Sequential()
model.add(LSTM(units=hp.Int('units1', min_value=32, max_value=256, step=32),
activation='tanh',
return_sequences=hp.Boolean('return_seq1', default=False),
input_shape=input_shape))
if hp.Boolean('dropout1', default=True):
model.add(Dropout(rate=hp.Float('dropout_rate1', 0.0, 0.5, step=0.1)))
if hp.Boolean('use_second_lstm', default=True):
model.add(LSTM(units=hp.Int('units2', min_value=16, max_value=128, step=16)))
if hp.Boolean('dropout2', default=True):
model.add(Dropout(rate=hp.Float('dropout_rate2', 0.0, 0.5, step=0.1)))
model.add(Dense(1))
lr = hp.Choice('lr', values=[1e-2, 1e-3, 1e-4])
optimizer = tf.keras.optimizers.Adam(learning_rate=lr, clipnorm=1.0)
model.compile(optimizer=optimizer, loss='mse')
return model
# --- Pipeline for LSTM Training and Prediction ---
def run_pipeline_lstm(data, test_ratio=0.2, lookback=10, max_trials=5, epochs=100):
X_train, y_train, dates_train, X_test, y_test, dates_test, scaler_X, scaler_y = prepare_data_lstm(data, target='Underlying_Close', test_ratio=test_ratio, lookback=lookback)
tuner = kt.RandomSearch(lambda hp: build_lstm_model(hp, input_shape=X_train.shape[1:]),
objective='val_loss',
max_trials=max_trials,
directory='lstm_tuner',
project_name='lstm_options_underlying')
tuner.search(X_train, y_train, epochs=epochs, validation_split=0.2,
callbacks=[tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=5)])
best_model = tuner.get_best_models(num_models=1)[0]
history = best_model.fit(X_train, y_train, epochs=epochs, validation_split=0.2,
callbacks=[tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=5)])
# Obtain predictions for both training and test sets
y_pred_train = best_model.predict(X_train)
y_pred_test = best_model.predict(X_test)
mse_train = mean_squared_error(y_train, y_pred_train)
mse_test = mean_squared_error(y_test, y_pred_test)
print("Train MSE:", mse_train)
print("Test MSE:", mse_test)
# Inverse transform the scaled targets
y_train_inv = scaler_y.inverse_transform(y_train)
y_pred_train_inv = scaler_y.inverse_transform(y_pred_train)
y_test_inv = scaler_y.inverse_transform(y_test)
y_pred_test_inv = scaler_y.inverse_transform(y_pred_test)
return best_model, scaler_X, scaler_y, tuner
if __name__ == '__main__':
features, options = build_black_scholes_dataset('MSFT', '2020-01-01', '2025-01-30')
print("Features:\n", features.tail())
print("\nOptions Data:\n", options.head() if not options.empty else "No options found")
best_model, scaler_X, scaler_y, tuner = run_pipeline_lstm(features, test_ratio=0.25, lookback=500, max_trials=10, epochs=200)